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O ■ ABSTRACT 

\ We present a measurement of the Lyman a flux probability distribution function 

(PDF) obtained from a set of eight high resolution quasar spectra with emission red- 
shifts at 3.3 < z < 3.8. We carefully study the effect of metal absorption lines on the 
shape of the PDF. Metals have a larger impact on the PDF measurements at lower 
Q ■ redshift, where there are relatively fewer Lyman a absorption lines. This may be ex- 

' plained by an increase in the number of metal lines which are blended with Lyman a 

c/2 . absorption lines toward higher redshift, but may also be due to the presence of fewer 

metals in the intergalactic medium at earlier times. We also provide a new measure- 
ment of the redshift evolution of the effective optical depth, r c g, at 2.8 < z < 3.6, and 

■ find no evidence for a deviation from a power law evolution in the log(r e ff ) — log(l + z) 

■ plane. The flux PDF measurements are furthermore of interest for studies of the ther- 
t-H \ mal state of the intergalactic medium (IGM) at z ~ 3. By comparing the PDF to 
£NJ . state-of-the-art cosmological hydrodynamical simulations, we place constraints on the 

■ temperature of the IGM and compare our results with previous measurements of the 
, PDF at lower redshift. At redshift z — 3, our new PDF measurements are consistent 
' with an isothermal temperature-density relation, T — TrjA 7-1 , with a temperature at 

O the mean density of T = 19250 ± 4800 K and a slope 7 = 0.90 ± 0.21 (la uncertain- 

£N| . ties). In comparison, joint constraints with existing lower redshift PDF measurements 

t-H ■ at z < 3 favour an inverted temperature-density relation with To = 17900 ± 3500 K 

^ \ and 7 = 0.70 ± 0.12, in broad agreement with previous analyses. 



Key words: cosmology: intergalactic medium - methods: numerical - quasars: ab- 
sorption lines. 



1 INTRODUCTION 

The Lyman a forest corresponds to the large number of ab- 
sorption features located blueward of the Lyman a emission 
line in the spectra of quasi-stellar objects (QSOs). In the 
last decade, thanks to the power of 10-m class telescopes 
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equipped with high-resolution spectrographs, considerable 
progress has been made toward understanding the nature of 
these absorption features. This population of discrete lines is 
due to HI absorption arising from the filamentary structure 
of the cosmic-web, and in the approximation of a Gaussian 
velocity dispersion can be well described by a series of Voigt 
profiles (Rauch 1998). A variable amount of absorption in 
the Lyman a forest is also due to ultraviolet (UV) transi- 
tions from heavy element ions which are usually associated 
with strong Lyman a lines. 

Quantities often studied in analyses of QSO spectra in- 
clude the power spectrum of the flux distribution, which is 
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used to constrain cosmological parameters and probe the 
dark matter power spectrum on scales of order of 50 h^ 1 
Mpc (e.g. Croft et al. 2002, Viel et al. 2004, McDonald et 
al. 2006). Another quantity thoroughly studied in the last 
few years is the flux probability distribution function (PDF), 
which is sensitive not only to the spatial distribution of dark 
matter, but also to the thermal state of the intergalactic 
medium (IGM). Previous studies of the Lyman a forest PDF 
have shed light on the physical state of the IGM mostly at 
z < 3. Bolton et al. (2008) compared the flux PDF measured 
from a sample of QSO spectra presented by Kim et al. (2007, 
hereafter K07) to a set of hydrodynamical simulations of 
the Lyman a forest, exploring different cosmological param- 
eters and various thermal histories. Agreement between the 
data and simulations was obtained by adopting a power-law 
temperature-density relation, T = ToA 7-1 (Hui & Gnedin 
1997), where the low density IGM (A = p/(p) < 10) was 
close to isothermal (i.e. with 7 ~ 1) or inverted (7 < 1). 
This result implied that the low-density regions of the IGM 
may be considerably hotter and their thermal state more 
complicated than usually assumed at these redshifts. 

Radiative transfer effects during the epoch of Hell 
reionization, which should occur at z > 3, are expected 
to play a non negligible role in setting this temperature- 
density relation. Indeed, evidence that the tail-end of Hell 
reionization is occurring at z ~ 3 is supported by a variety 
of observations, such as the opacity of the Hell Lyman a for- 
est (Reimers et al. 1997; Heap et al. 2000; Kriss et al. 2001, 
Shull et al. 2010, Worseck et al. 2011, Syphers et al. 2011), 
an increase in the temperature of the IGM at mean density 
toward z ~ 3 (Schaye et al. 2000, Theuns et al. 2002, Lidz 
et al. 2010, Becker et al. 2011, but see McDonald et al. 2001, 
Zaldarriaga, Hui, & Tegmark 2001), the observed evolution 
of the CIV/SilV metal-line ratio (Songaila & Cowie 1996; 
Songaila 1998, but see Kim, Cristiani & D'Odorico 2002) 
and the evolution of the Lyman a forest effective optical 
depth (Bernardi et al. 2003, Faucher-Giguere et al. 2008, 
Paris et al. 2011, but see Bolton, Oh & Furlanetto 2009b). 

However, recent work suggests that photo-heating dur- 
ing Hell reionization alone is not sufficient to achieve an in- 
verted temperature-density relation, 7 < 1 (McQuinn et al. 
2009; Bolton, Oh & Furlanetto 2009a). A recently proposed 
alternative channel for heating, which naturally produces 
an inverted temperature-density relation, may be provided 
by very energetic gamma rays from blazars (Chang et al 
2011, Puchwein et al. 2011). In this scenario, very high en- 
ergy gamma rays annihilate and pair produce on the extra- 
galactic background light. Powerful plasma instabilities then 
dissipate the kinetic energy of the electron-positron pairs lo- 
cally, volumetrically heating the IGM. On the other hand, 
continuum placement errors could partly explain the prefer- 
ence of the K07 data for an inverted temperature-density re- 
lation (see Lee 2012 for a recent discussion). Indeed, Bolton 
et al. (2008) and more recently Viel et al. (2009) found 
an isothermal temperature-density relation was consistent 
(within 1.5a) with the K07 PDF constraints when treating 
continuum uncertainties conservatively. 

Ultimately, however, for further progress in understand- 
ing the origin of an inverted temperature-density relation, 
additional observational constraints on the IGM thermal 
state are required. An investigation of the flux probabil- 
ity distribution at redshifts around z ~ 3 or higher may 



therefore shed more light on the possible inversion of the 
temperature-density relation. Such a measurement is also 
complementary to other constraints on the IGM thermal 
state, such as those presented recently by Becker et al. 
(2011), where the temperature-density relation is not mea- 
sured directly. Such an investigation is the main motivation 
of the present work. 

In this paper, we use eight high resolution QSO spectra 
observed at 3.3 < z < 3.8 to extend the study of the Ly- 
man a forest PDF by K07 to higher redshifts. The aim of 
this paper is two-fold. First, we perform a detailed study of 
the redshift evolution of the flux PDF across a wide redshift 
range. This enables us to study possible systematic trends 
in the PDF due to redshift-dependent effects, such as the 
presence of metal absorption lines. In order to corroborate 
our analysis, we also measure the Lyman a effective optical 
depth and compare it to different estimates by other authors. 
In the second step, the measured flux PDF is compared to 
theoretical results obtained from cosmological hydrodynam- 
ical simulations. This enables us to investigate the thermal 
state of the IGM during this interesting redshift range. 

This paper is organized as follows. In Section 2, we 
present the QSO spectra and the data analysis, and in Sec- 
tion 3 we present our measurements of the flux PDF and 
effective optical depth. In Section 4, we introduce the hy- 
drodynamical simulations used to obtain constraints on the 
thermal state of the IGM from the PDF and, in Section 5, 
we report the results from the analysis of the simulations. 
Finally, in Section 6, we present our conclusions. 



2 THE DATA SET 

The data set used in this paper consists of eight high- 
redshift QSO spectra. Two of the spectra (PKS2126-158 
and PKS2000-330) were obtained with the Ultraviolet and 
Visual Echelle Spectrograph (UVES) (Dekker et al. 2000) 
at the Kueyen unit of the European Southern Observa- 
tory (ESO) VLT (Cerro Paranal, Chile) as part of the ESO 
Large Programme (LP): 'The Cosmic Evolution of the IGM' 
(Bergeron et al. 2004). PKS1937-101 is from the program 
Carswell et al. 077.A-0166(A), Q1317-0507, SDSSJ16212- 
0042 and Q1249-0159 are from Kim et al. 075.A-0464(A), 
Q1209+0919 is from the program Dessauges-Zavadsky et al. 
073.B-0787(A) and Q0055-269 is from D'Odorico et al. 65.0- 
0296(A). All the spectra have been reduced with the stan- 
dard UVES pipeline provided by ESO. Wavelengths have 
been corrected to vacuum-heliocentric. 

The properties of the QSO spectra used in this paper 
are summarized in Table 1. The velocity resolution of each 
spectrum is 6.7 km s _1 full-width half maximum, and each 
spectrum is binned in 0.05 A pixels. The signal-to-noise 
(S/N) ratio has been calculated by considering regions of 
the spectra redwards of the QSO emission peak, and corre- 
sponds to the inverse of the standard deviation of the aver- 
age normalized flux for pixels free of absorption. For some 
QSOs, the spectra have already been used by D'Odorico et 
al. (2010) for the study of the evolution of CIV in the IGM. 
The spectra of three QSOs (Q1317-0507, Q1249-0159 and 
SDSSJ16212-0042) have been reduced for this work. 

Most of the QSOs include one or more Lyman Limit 
systems (LLS), i.e. absorption systems with neutral col- 
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Table 1. The QSO sample used in this work. In the first column the name of the QSO is reported. Subsequent columns list the emission 
redshift of the QSO, the redshift range covered by the Lyman a forest, the wavelength range covered by the Lyman a forest and the 
signal to noise (S/N, see text for details). Note the wavelength range of the Lyman a forest excludes the region 4000 km s _1 blucwards 
of the Lyman a emission line in order to avoid the proximity effect (see Sect. 2. 2). 



QSO 


^ a 
~em 




A Lya (A) 


S/N 


Notes 






PKS2126-158 


3.292 


2.68-3.23 


4473-5148 


110 


1 


LLS 




=2.967), 












2 


Sub-DLAs (z=2.638, 2.769) 


Q1209+0919 


3.295 


2.62-3.23 


4404-5150 


25 


8 


LLS 






Q0055-269 


3.66 


2.93-3.59 


4775-5583 


51 










Q1249-0159 


3.66 


2.94-3.60 


4784-5594 


70 


3 


LLS 


(* 


=3.102, 3.5243, 3.5428) 


SDSSJ16212-0042 


3.70 


2.97-3.63 


4821-5637 


93 


2 


LLS 


(*= 


=3.106, 3.144) 


Q1317-0507 


3.71 


2.98-3.66 


4840-5659 


86 


1 


LLS 


(*= 


=3.288) 


PKS2000-330 


3.78 


3.03-3.72 


4903-5733 


68 


2 


LLS 


(z-- 


=3.189, 3.33) 


PKS1937-101 


3.79 


3.04-3.72 


4910-5742 


150 











° In some cases, the emission redshift z cm estimated on the basis of the QSO Lyman a emission line is an underestimation of the true 
value (see K07). If Lyman a absorption lines with z > z em are present, we assume the highest redshift of the absorption line as a proxy 
for z em . 



1 

0.5 

1 

0.5 


1.5 
1 

0.5 


1.5 
1 

5 


1 

5 

1 

0.5 

1 

0.5 

1 

0.5 




PKS2 




58 - mask 



Q1209+0919 -lobs 




: Q0055 269 obs 

iiflfiiiiif 



Q0055-269 - mask 



Q 1249-0 159 - obs. 

mmmnmm 





Q1249 0159 mask 



SDSSJ16212-0042 - obs 





S )SSJ 16212-0042 - mask 



TT1I 



i i i i { r i 

Q1317-0507 - obs. 



Q1317-0507 - mask 



:'KS 2000 130 obs. 



1000 1050 1100 1150 1200 

A(A) - Rest frame 



1050 1100 1150 

X(A) - Rest frame 




1 

0.5 



1 

0.5 



1 

0.5 



1 

0.5 




: 200 



Figure 1. The eight QSO spectra which comprise our sample. Each spectrum is normalised to the quasar continuum and the wavelengths 
are displayed in the QSO rest-frame. For each QSO, the upper panel shows the normalised spectrum, whereas the lower panel displays 
all the spectral regions containing metal lines. To measure the PDF, these regions have been masked following the criteria described in 
Sect. 2.3. In each panel, the left dotted line and the right dotted line indicate the wavelengths of the Lyman /3 and Lyman a transitions. 
The regions with no flux, such as the one apparent in the spectrum of PKS 2000-330 centered at A ~ 1180 A, are due to the wavelength 
gaps caused by using three separate CCDs in the dichroic setting. 
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umn density 17.2 < log^m/cm -2 ) < 19. For the com- 
parison of our data to simulation results, we have therefore 
excised all the pixels belonging to the spectral regions in- 
cluded in the range Alls — 50 A < A < Alls + 50 A, 
where Alls is the central wavelength of the LLS. PKS2126- 
158 additionally includes two sub-damped Lyman a systems 
(19 < log(7V H i/cm- 2 ) < 20.3). The spectral regions where 
these systems are present have also been excised from the 
spectrum, following the same criterion as for LLS. In addi- 
tion, the identified metal lines associated with both LLS and 
sub-DLAs have been removed from the spectra as described 
in Sect. 2.3. 

The two QSOs PKS2126158 and Q0055-269 were also 
included in the K07 sample. It is worth noting that the wave- 
length ranges for PKS2126158 and Q0055-269 used in this 
work are different to those reported in Tab. 1 of K07. In the 
case of PKS2126158, our blue wavelength cut is different 
to the K07 value due to the presence of two sub-DLAs at 
z — 2.638 and z = 2.769. We excluded the two regions con- 
taminated by the sub-DLAs as described above, and a small 
region between the two sub-DLAs was left in our spectrum, 
extending from 4470 A to 4540 A. K07 instead decided to re- 
move this region and cut their spectrum at 4630 A. The red 
wavelength cut was chosen by K07 in order to avoid having 
pixels at wavelengths > 5112 A corresponding to redshifts 
z > 3.2. However, this redshift range is important for study- 
ing the possible effects of Hell reionization and so we include 
it in this work. In the case of Q0055-269, the blue cut was 
performed in our spectrum at a wavelength corresponding 
to the Lyman f} emission line, whereas in the case of K07 
they performed it at 4785 A in order to exclude a metal 
absorption system at 4784 A. The red cut was performed by 
K07 at 5112 A to again discard pixels with z > 3.2. 

In Fig. 1, we show the reduced spectra for all the QSOs 
of our sample and, for each spectrum, the contamination 
from metal lines. The regions containing metal lines have 
been masked following the criteria described in Sect. 2.3. In 
Fig. 2, for each QSO spectrum, the redshift coverage and 
the Lyman a forest wavelength coverage are shown. 

2.1 Continuum fitting 

In each spectrum all the absorption lines redwards of the 
Lyman a emission line, as well as those belonging to the 
Lyman a forest, have been fitted with Voigt profiles. To 
satisfactorily perform this task, for each spectrum a precise 
estimate of the continuum is required. However, the problem 
of continuum fitting high resolution data requires particu- 
lar attention, and it has been widely discussed in the lit- 
erature (e.g. K07; Desjacques et al. 2007; Faucher-Giguere 
et al. 2008). An automatic procedure to fit QSO spectra 
has also been developed by Dall' Aglio et al. (2008a; 2008b) . 
They used an extended sample of high resolution, high S/N 
QSO spectra, and by means of an adaptive techinque involv- 
ing cubic spline interpolations and a Monte-Carlo method, 
they were able to estimate the QSO continua and relative 
uncertainties. In this paper, however, continuum fitting was 
performed on our smaller data set as described below. It will 
be interesting to use a larger sample in the future to develop 
automatic techniques for continuum fitting and investigate 
associated uncertainties further. 

For each QSO, the continuum was fitted using the IRAF 
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Figure 2. The solid lines display the redshift (and wavelength, 
see upper axis) coverage of the Lyman a forest in the QSO spectra 
used in our sample. The dashed lines at z = 2.62, z = 3.17 and 
z = 3.72 represent the extremes of the two redshift bins used to 
calculate the flux probability distribution function, as described 
in Sect. 3. 



package, by means of high-order polynomial functions. The 
continuum fitting was performed by cutting each spectrum 
into several chunks, depending on the S/N and on the num- 
ber of absorption lines, and then joining the chunks back 
together into one single spectrum. The continua were de- 
termined by using regions free of absorption lines. This task 
does not present substantial difficulties when the region red- 
wards of the QSO Lyman a emission line is considered, 
since this region presents only a few (metal) absorption lines. 
However, regions bluewards of the QSO Lyman a emission 
lines, such as the Lyman a forest, are instead populated by 
a large number of absorption lines. In general, the higher the 
emission redshift of the QSO, the greater the average Lyman 
q forest opacity. This makes continuum fitting toward higher 
redshift more difficult, since regions which are free of absorp- 
tion lines become increasingly small and rare. In this case, a 
linear interpolation between two subsequent absorption- free 
chunks is used as a first-order guess for the continuum level. 
Another major problem in continuum fitting concerns the 
possible presence of wide and shallow absorption features in 
the spectra. Large-scale absorption can sometimes be incor- 
rectly attributed to unabsorbed parts of the spectrum. This 
represents a problem in this procedure, since it can lead to 
an underestimate of the true continuum. To avoid such mis- 
takes, it is necessary not only to check carefully every small 
chunk of the spectrum which is being fitted, but also the 
joint spectrum by searching for possible wide depressions 
due to large-scale absorption features which are not visible 
while working on the individual chunks. All the spectra and 
the fitted continua have been repeatedly and independently 
checked in this manner in order to identify such features and 
to avoid continuum underestimates. 

As discussed in K07, the uncertainty on the continuum 
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Figure 3. An illustration of the uncertainty on the continuum placement for three of the spectra in our sample with different S/N 
values: S/N=25 for Q1209+0919 (upper panel), S/N=68 for PKS2000-330 (middle panel) and S/N=150 for PKS1937-101 (lower panel). 
In each panel, the final continuum fit is shown by the solid curve, while the dotted and dashed curves show the fitted continuum after it 
is increased/decreased by 5 per cent and by 10 per cent, respectively. 



depends on the S/N, which varies not only from spectrum 
to spectrum, but may even vary within one single spec- 
trum. As a general rule, the higher the S/N, the smaller the 

continuum uncertainty. This is illustrated in Fig. 3^ where 
we show three unnormalized spectra with their final fitted 
continuum (solid lines), and with the continuum increased 
and decreased by 5 per cent (dotted lines) and by 10 per 
cent (dashed lines). The cases depicted in Fig. 3 represent 



t In Fig. 3, the spectrum of PKS1937-101 appears to be nois- 
ier than PKS2000-330, contrary to what is reported in Tab. 1. 
Recall, however, that the S/N reported in Tab. 1 have been 
calculated by means of absorption-free regions redwards of the 
QSO emission peak, and that S/N may increase with increas- 
ing wavelengths owing to the increase in the efficiency of the 
UVES CCDs. As a check, we have therefore recalculated the S/N 
in the Lyman a forest for PKS1937-101 and PKS2000-330 by 
computing the average ratio between the flux and relative er- 
ror < F/ap >, finding in both cases < F/ap >= 59. How- 
ever, the maximum S/N value found in the Lyman a forest of 
PKS1937-101 is < F/ap > ma x^ 100, whereas for PKS2000-330 
it is < F/ap > ma x— 62. 



the spectrum with the lowest S/N in Tab. 1 (Q1209+0919, 
upper panel), one case with an intermediate S/N value 
(PKS2000-330, middle panel) and the spectrum with the 
highest S/N (PKS1937-101, lower panel). Although it is very 
difficult to precisely quantify the uncertainty for the contin- 
uum determination, Fig. 3 demonstrates that, in the case of 
Q1209+0919, an uncertainty of the order of 10 per cent is 
exaggerated, but an uncertainty of 5 per cent seems plau- 
sible. As stressed earlier, this case represents the spectrum 
with the worst S/N in our sample. In the other two cases, 
it appears that a continuum which is decreased or increased 
by 5 per cent represents an underestimation and overestima- 
tion, respectively, of the actual continuum. For this reason, 
a typical uncertainty of 5 per cent should be regarded as 
conservative and representative of most of the spectra in 
our sample, in agreement with the analysis of K07 and the 
one of DalPAglio et al. (2008a). Later on, the effects of the 
continuum estimate on the flux PDF will be investigated in 
more detail. 
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2.2 Voigt Profile fitting 

Once the continuum fitting was completed, in each spec- 
trum, all the absorption features blueward of the QSO Ly- 
man q emission line were then fitted with Voigt profiles. The 
aim of this procedure is to determine, for each absorption 
line, the redshift z, the Doppler parameter b and the col- 
umn density N of the physical system giving rise to that 
absorption feature. 

The Lyman a forest includes the spectral region be- 
tween the wavelengths of the Lyman /3 emission line, Ap = 
1025.72 (1 + zqso), and the wavelength corresponding to the 
Lyman a emission line, i.e. X a — 1215.67(1 + zqso)- To 
avoid contamination from the QSO line-of-sight proximity 
effect, in each spectrum we excluded an interval of 4000 km 
s~ x bluewards of the Lyman a emission line. All the absorp- 
tion features were then identified and fitted by means of the 
RDGEN§, VPFIT^ and VPGUESsll packages. 

Firstly, all the metal absorption lines redwards of the 
QSO Lyman a emission peak were identified and fitted. This 
task does not present significant difficulties, since the num- 
ber of absorption lines in this spectral region is relatively 
low and lines due to different transitions are rarely blended. 
Secondly, the identification of heavy element absorption in 
the Lyman a forest was performed. This task requires more 
attention owing to the large amount Lyman a lines in this 
region and to the frequent presence of broad, saturated ab- 
sorption features. In this spectral region, isolated metal lines 
are very rare. Most of the metal absorption lines are blended 
with Lyman a lines, complicating the precise identification 
of metal line profiles. A common technique to search for 
metal absorption in the Lyman a forest therefore employs 
the use of velocity width profiles. Once a strong transition 
is identified at wavelengths redward of the QSO Lyman a 
emission line, e.g. due to CIV or SilV, associated ionic tran- 
sitions such as Sill, Silll, CII at the corresponding redshift 
and velocity profile can be searched for in the Lyman a for- 
est. 

When metal lines in the Lyman a forest are heavily 
blended with HI absorption, it is generally difficult to iden- 
tify the shape of the profile. The parameters which are most 
difficult to recover in this case are the Doppler parameter 
b and the column density N. In these cases, as an initial 
guess for VPFIT we provide the profile of the lines asso- 
ciated with the systems already identified redwards of the 
QSO emission, usually satisfactorily fitted by VPFIT. VP- 
FIT then corrects the initial guess and modifies it in order 
to obtain an improved fit to that spectral region. Once a sta- 
ble solution is obtained which is not sensitive to the initial 
guess and characterised by a satisfactory reduced x 2 value 
(typically of the order of 1 — 3, depending on the noise level 
of the spectrum) the fit process is complete. An alterna- 
tive way to identify metal absorption in the forest is based 
on the value of the Doppler parameter, which is in general 
narrower than for Lyman a lines (e.g. Tescari et al. 2011). 
However, in this analysis any absorption not corresponding 



§ http://www.ast.cam.ac.uk/ rfc/rdgen.html 
^ http://www.ast.cam.ac.uk/ rfc/vpfit.html 
II http://www.eso.org/ jliske/vpguess/ 
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Figure 4. The flux PDF errors as a function of the normalised 
flux measured from all eight QSOs, computed by means of the 
jackknife (black solid lines) and the bootstrap (red dashed line) 
methods. 

to an identified heavy element transition has been assumed 
to be due to Lyman a absorption. 

2.3 Metal removal from the spectra 

The presence of metal absorption alters the transmitted frac- 
tion in the QSO spectra and consequently the shape of the 
PDF. For this reason, in order to use the PDF as a probe 
of the physical state of the Lyman a forest, it is necessary 
to remove the metal absorption from the Lyman a forest in 
the QSO spectra. 

Once the metal identification and Voigt profile fitting 
is complete, for each QSO we remove metal lines from the 
spectra in the following manner. The pixels belonging to 
spectral regions contaminated by identified metal absorp- 
tion are excised from the spectra using a quantitative crite- 
rion which depends on the Doppler parameter of the metal 
absorption line, b. All the pixels at wavelengths, A, within 
the interval Ao — AA < A < Ao + AA were removed from 
the spectra, where Ao is the central wavelength of the metal 
absorption line and AA = Xz,o (1 + 2 ')- Here Xz,o is the rest- 
frame wavelength of the transition of a given element Z and 
z' — z — dz, where z is the redshift of the absorbing system 
and dz = 2b ln ^( 1 + z \ p or eac h QSO spectra, we created a 
"mask" using this procedure, which is then used to remove 
the pixels belonging to the metal-contaminated regions. 

The completeness of the masking is dependent on the 
success of identifying metals. It is clearly impossible to iden- 
tify all the metal absorption in the Lyman a forest. In this 
work we are dealing with QSO spectra at z > 3 which con- 
tain a large number of HI absorption lines. At these red- 
shifts, metal lines are frequently blended with strongly sat- 
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urated HI absorption, and when this is the case their contri- 
bution to the total absorption is negligible. However, there 
will be very few cases where the metal lines are isolated and 
not identified. For this reason, we believe that the results 
presented in this paper are not substantially affected by in- 
complete metal identification in the Lyman a forest. This 
issue will be discussed further in the section dedicated to 
the results. 



3 DATA ANALYSIS AND RESULTS 

Once each spectrum has been normalised to the continuum, 
i.e. once the quantity F — S/C has been calculated for each 
pixel, where S and C are the observed flux and the esti- 
mated continuum, respectively, the flux probability distribu- 
tion function can be measured. For each flux bin F, the PDF 
is calculated from the ratio between the number, JVf,f+af, 
of pixels with flux between F and F + AF and the total 
number, AT tot , of pixels in the spectral region corresponding 
to the Lyman a forest (K07): 

PDF(F) = ^±^. (1) 

JVtot 

The PDF is measured using this procedure in bins of width 
AF = 0.05 (McDonald et al. 2000; K07). We have assigned 
flux values F = and F = 1 to all the pixels characterised 
by flux levels F < 0.025 and F > 1, respectively. The PDF 
is calculated both for single QSO spectra, and for the full 
sample. To study the redshift evolution of the PDF, we also 
divide the full sample in two redshift bins. The redshift bin 
width is chosen in order to have comparable numbers of 
pixels in each bin. 



3.1 Error bar estimates 

The error bars are estimated by means of a jackknifc 
method, as described by Lidz et al. (2006). Once the PDF 
for a single QSO or for the full sample in a given redshift 
bin has been calculated, the contributing spectra are then 
divided into n c chunks of width ~ 50 A. If the PDF in the 
flux bin F\ is P(Fi) and the PDF estimated without the k- 
th chunk in the flux bin Fi is Fk(F), the covariance matrix 
cov(i,j) between the PDF in flux bin Fi and the PDF in 
flux bin Fj can be defined as 

cov(i, j) = ^[P(F ; ) - P k (F i )][P(F j ) - P k (Fj)], (2) 

k=l 

For the flux bin Fi, the errors computed with the jackknifc 
method, <t.j, are given by the square root of the diagonal of 
the matrix cov(i,j), i.e. 

a j = >/cov(i,i). (3) 

In addition, it is useful to have an alternative way to es- 
timate the error bars for comparison to the uncertainties 
calculated using the jackknife method. We therefore also es- 
timate error bars with a bootstrap technique. This technique 
consists of dividing the spectra into chunks of ~ 5 A , and 
drawing a large number (500) of these chunks from the sam- 
ple randomly with replacement. The PDF averaged over the 



500 realizations coincides with the PDF of the random sam- 
ple. For each flux bin, the uncertainty is then estimated by 
calculating the standard deviation of the PDFs obtained by 
sampling the 500 chunks with replacement. 

In figure 4, the la errors for the PDF, estimated using 
the jackknife and bootstrap techniques, are compared. The 
PDF is computed from the full sample in this instance, with- 
out dividing the spectra into different redshift bins. Fig. 4 
shows how the two methods provide comparable error esti- 
mates, although the jackknife estimates are larger at F = 
and F = 1. In the remainder of the paper, unless otherwise 
stated, the uncertainties for the PDF will be computed by 
using the jackknife method. 

3.2 The effect of the metals 

The effect of metal absorption on the PDF is displayed in 
Figures 5 and 6, where we show the PDF obtained from 
each of the eight QSOs in our sample. In these two figures, 
the solid and dashed curves in the large panels show the 
PDF computed including and excluding metal absorption 
lines, respectively. In the small panels, for each QSO the 
ratio between the PDF calculated including metals and after 
removing the metals is shown. The average Lyman a forest 
redshift spans from 2.96, for PKS2126-158, up to 3.38 in the 
case of PKS1937-101. 

In general, the effect of removing metals on the PDF is 
to decrease the number of pixels with low flux and increase 
the number of pixels with high flux. As a result, the PDF 
computed after metal removal is slightly steeper than the 
PDF calculated before removing the metal absorption. This 
effect is most clearly visible for some of the QSOs in Fig. 5, 
such as PKS2126-158. For this spectrum, the effect of metal 
removal on the PDF is visible for flux values F < 0.6, and 
is in agreement with the lower redshift analysis of K07, who 
find a PDF ratio similar to ours. 

The effect of metal absorption on the PDF is also ex- 
pected to depend on the redshift of the QSO and hence the 
average redshift of the Lyman a absorption. The impact of 
metal absorption on the PDF is expected to be more sig- 
nificant at lower redshifts, where Lyman a forest opacity is 
generally lower (K07). A weak trend with redshift is indeed 
visible in Figs. 5-6, but it is illustrated more clearly in Fig. 7, 
where for each QSO we show the ratio between the PDF cal- 
culated including metal lines and the PDF computed after 
removing metals from the spectra. The ratio is higher than 
1.2 for flux values F < 0.4 and lower than 0.9 for F > 0.9 
only in one case, PKS2126-158, which is one of the QSOs 
with the lowest value of (z) . Q1209+0919 has a compa- 
rable (z) and has a ratio close to 1.2 at F < 0.1. The other 
six QSOs have similar (z) (between 3.26 and 3.38), and also 
exhibit comparable levels of metal contamination. For these 
six QSOs small fluctuations from one LOS to another are 
apparent, but the ratios are all between 0.9 and 1.2 and are 
much smaller than the variations for the QSOs at (z) < 3. 

The redshift dependence of metal absorption may be 
understood in terms of the larger number of Lyman a ab- 
sorption lines toward higher redshift, which leave less room 

** For each QSO, (z) is the average redshift of the Lyman a 
forest range as reported in Tab. 1. 
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Figure 5. The effect of metal line absorption on the PDFs obtained for four of the QSOs in our sample. In the main panels, the solid 
curves display the PDFs calculated for each spectrum before removing metal lines. The solid squares with error bars show the PDFs 
calculated for the masked spectra, i.e. after removing the metal absorption. The dashed lines are drawn to guide the eye through the 
solid squares. For each QSO, the name and the average redshift of Lyman a forest absorption is indicated. In the smaller panels, the 
solid curve represents the ratio between the PDF calculated before and after removing metal lines. 
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Figure 6. As in Fig. 5, but now for the other four QSOs in our sample. 
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Table 2. The mean PDF measured from the full sample, divided into two redshift bins and calculated by (1) using all pixels, (2) after 
removing contamination from metal absorption lines and (3) after removing Lyman limit systems and metal lines. Note that the redshift 
bins in the final two columns have a slightly different mean redshift due to the removal of LLS. The redshift bins have been slightly 
modified in order to have comparable numbers of pixels in each bin. 





All pixels 


No metals 


No metals, No LLS 




(1) 


(2) 


(3) 


F 


< 2 > = 2.90 


< z > = 3.45 < z > = 2.90 


<z> = 3.45 <z>=2.95 < z > = 3.48 



0.00 


1.76 


± 


0.18 


2.37 


± 


0.24 


2.00 


± 


0.16 


2.38 


± 


0.23 


1.54 


± 


0.11 


2.28 


± 


0.30 


0.05 


0.66 


± 


0.05 


0.71 


± 


0.04 


0.74 


± 


0.04 


0.72 


± 


0.04 


0.57 


± 


0.04 


0.70 


± 


0.06 


0.10 


0.41 


± 


0.03 


0.41 


± 


0.03 


0.46 


± 


0.03 


0.43 


± 


0.03 


0.34 


± 


0.03 


0.42 


± 


0.04 


0.15 


0.36 


± 


0.02 


0.34 


± 


0.02 


0.39 


± 


0.02 


0.36 


± 


0.02 


0.30 


± 


0.03 


0.35 


± 


0.03 


0.20 


0.30 


± 


0.01 


0.33 


± 


0.02 


0.32 


± 


0.01 


0.35 


± 


0.02 


0.26 


± 


0.02 


0.35 


± 


0.02 


0.25 


0.30 


± 


0.02 


0.33 


± 


0.02 


0.32 


± 


0.02 


0.34 


± 


0.02 


0.26 


± 


0.03 


0.35 


± 


0.03 


0.30 


0.30 


± 


0.02 


0.36 


± 


0.02 


0.32 


± 


0.02 


0.37 


± 


0.02 


0.27 


± 


0.02 


0.36 


± 


0.03 


0.35 


0.33 


± 


0.02 


0.35 


± 


0.02 


0.35 


± 


0.02 


0.36 


± 


0.02 


0.30 


± 


0.02 


0.35 


± 


0.02 


0.40 


0.30 


± 


0.02 


0.37 


± 


0.02 


0.34 


± 


0.02 


0.38 


± 


0.02 


0.29 


± 


0.03 


0.36 


± 


0.03 


0.45 


0.31 


± 


0.01 


0.41 


± 


0.01 


0.34 


± 


0.01 


0.42 


± 


0.02 


0.31 


± 


0.02 


0.42 


± 


0.02 


0.50 


0.35 


± 


0.02 


0.46 


± 


0.01 


0.38 


± 


0.02 


0.47 


± 


0.02 


0.34 


± 


0.02 


0.48 


± 


0.02 


0.55 


0.40 


± 


0.02 


0.58 


± 


0.02 


0.43 


± 


0.02 


0.58 


± 


0.02 


0.40 


± 


0.04 


0.60 


± 


0.03 


0.60 


0.51 


± 


0.02 


0.61 


± 


0.03 


0.52 


± 


0.02 


0.62 


± 


0.03 


0.50 


± 


0.03 


0.62 


± 


0.04 


0.65 


0.56 


± 


0.02 


0.71 


± 


0.03 


0.58 


± 


0.02 


0.72 


± 


0.04 


0.57 


± 


0.03 


0.72 


± 


0.05 


0.70 


0.70 


± 


0.04 


0.81 


± 


0.02 


0.72 


± 


0.04 


0.82 


± 


0.03 


0.68 


± 


0.03 


0.84 


± 


0.05 


0.75 


0.91 


± 


0.05 


0.95 


± 


0.03 


0.91 


± 


0.05 


0.96 


± 


0.04 


0.86 


± 


0.04 


0.97 


± 


0.04 


0.80 


1.21 


± 


0.05 


1.13 


± 


0.06 


1.19 


± 


0.05 


1.14 


± 


0.05 


1.14 


± 


0.05 


1.15 


± 


0.06 


0.85 


1.55 


± 


0.06 


1.37 


± 


0.06 


1.49 


± 


0.06 


1.38 


± 


0.05 


1.44 


± 


0.07 


1.40 


± 


0.06 


0.90 


2.01 


± 


0.09 


1.88 


± 


0.08 


1.92 


± 


0.08 


1.87 


± 


0.08 


2.01 


± 


0.10 


1.93 


± 


0.10 


0.95 


2.84 


± 


0.07 


2.66 


± 


0.13 


2.67 


± 


0.07 


2.56 


± 


0.14 


3.24 


± 


0.12 


2.69 


± 


0.21 


1.00 


3.93 


± 


0.19 


2.89 


± 


0.26 


3.61 


± 


0.18 


2.77 


± 


0.27 


4.38 


± 


0.29 


2.69 


± 


0.24 



Table 3. The HI effective optical depth and its dependence on the contamination of metal absorption lines, calculated after dividing 
each spectrum in two redshift bins. In columns 1, 2 and 3 the name of the QSO, the Lyman a wavelength range and the average redshift 
are shown, respectively. In columns 4, 5 and 6 we report the HI effective optical depth with its error, the optical depth computed without 
removing metal lines with its error, and the percentage contribution of metal absorption, respectively. 



QSO wavelengths (A) <z> r off T cffjH i +z "r cffiHI+z /"r off metals (per cent) 



PKS2126-158 


4470-4880 


2. 


85 





329 


± 


0. 


037 





.421 


± 





.041 


1 


.278 


± 


0. 


190 


28 




4880-5148 


3 


.12 





.238 


± 





.031 





.271 


± 





.031 


1 


.138 


± 


0. 


197 


14 


Q1209+0919 


4404-4780 


2. 


78 





373 


± 


0. 


031 





.423 


± 





.037 


1 


.132 


± 


0. 


136 


13 




4780-5150 


3. 


.08 





.425 


± 





.040 





.450 


± 





.039 


1 


.060 


± 


0. 


135 


6 


Q0055-269 


4774-5170 


3 


09 





339 


± 


0. 


033 





.344 


± 





.032 


1 


.014 


± 


0. 


137 


1 




5170-5583 


3 


.42 





.429 


± 





.036 





.439 


± 





.035 


1 


.024 


± 


0. 


119 


2 


Q1249-0159 


4784-5190 


3 


10 





382 


± 


0. 


028 





.399 


± 





.028 


1 


.045 


± 


0. 


106 


4 




5190-5594 


3 


44 





.529 


± 





.041 





.542 


± 





.040 


1 


.024 


± 


0. 


110 


2 


SDSSJ16212-0042 


4821-5235 


3. 


14 





555 


± 


0. 


042 





.596 


± 





044 


1 


.073 


± 


0. 


113 


7 




5230-5637 


3 


.47 





.394 


± 


0. 


.032 





.427 


± 





.034 


1 


.085 


± 


0. 


123 


8 


Q1317-0507 


4840-5240 


3. 


15 


0. 


.384 


± 


0. 


031 





.388 


± 





031 


1 


.010 


± 


0. 


116 


1 




5240-5659 


3. 


.48 





.495 


± 





.038 





.499 


± 





.039 


1 


.008 


± 


0. 


110 


1 


PKS2000-330 


4903-5270 


3. 


18 





574 


± 


0. 


048 





.582 


± 


0. 


.046 


1 


.015 


± 


0. 


117 


1 




5270-5734 


3 


.53 





.553 


± 





.043 





.570 


± 





.041 


1 


.029 


± 


0. 


108 


3 


PKS1937-101 


4910-5323 


3. 


21 





.371 


± 


0. 


027 





.376 


± 





.027 


1 


.014 


± 


0. 


103 


1 




5323-5742 


3 


55 





.487 


± 





.037 





.492 


± 





.036 


1 


.011 


± 


0. 


106 


1 



for isolated metal absorption features. This does not neces- 
sarily mean that QSOs at higher redshifts have less metal 
absorption systems along their LOS, as the metal lines may 
instead be blended with Lyman a absorption lines. At the 
same time, however, this observation does not fully rule out 
the possibility the metal content of the intergalactic medium 



is decreasing toward higher redshift. We will return to the 
issue of metal contamination of QSO spectra in Sect. 3.5. 
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Figure 7. The ratios between the PDFs measured for each of 
the eight QSO spectra analysed in this work before and after the 
removal of metal line absorption. Note that the line thickness 
increases with the average Lyman a forest redshift. 



3.3 The effect of the continuum uncertainty 

The effects of the systematic uncertainty due to the con- 
tinuum fitting procedure are shown in Fig. 8 for the full 
sample in the redshift bin centered at z = 2.9. In this fig- 
ure, the PDF is measured without removing the metal lines 
and is shown for the final continuum level (solid curve), as 
well as the PDF obtained by increasing and decreasing the 
continuum level by 5 per cent (dashed and dotted curves, 
respectively). Here the continuum is conservatively rescaled 
by 5 per cent since we have seen from Fig. 3 this value is 
representative of the continuum uncertainty for most of the 
spectra in our sample. In both cases, the strongest effects on 
the PDF are visible at high flux values (F > 0.5), whereas 
a rescaling of the continuum has a negligible effect in the 
lowest flux bins. 

The net effect of increasing the continuum by 5 per cent 
is to remove pixels at F ~ 1, moving the peak of the original 
PDF to lower flux bins. The bins at values F < 0.9 are also 
populated with more pixels than is the case for the final 
continuum, and the slope of the PDF between F ~ 0.1 and 
F ~ 0.9 tends to increase overall. On the other hand, by 
decreasing the continuum by 5 per cent, pixels are instead 
moved towards higher flux bins and the PDF peaks at F = 1. 

The determination of the continuum is a very prob- 
lematic task and subject to several uncertainties. However, 
Fig. 8 shows how the PDF would vary if we had systemat- 
ically underestimated or overestimated the continuum by 5 
per cent for all absorption systems. This seems rather im- 
plausible: it is far more likely that our continuum fitting 
method underestimates or overestimates the true continuum 
in chunks of our spectra, rather than consistently across the 
entire data set.tt If this is the case, the effect of the contin- 
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Figure 8. The effect of continuum placement uncertainty on the 
PDF calculated for our full sample in the lower redshift bin. In 
the upper panel, the solid line is the PDF calculated with the 
final continuum. The solid squares with error bars are the PDF 
values in each flux bin with the associated error estimated by 
means of the jackknife method described in Sect. 3.1. The short- 
dashed line, the long-dashed line and the dot-long-dashed line 
represent the PDF calculated by increasing the continuum level 
by 5, 3 and 1.5 per cent, respectively. The dotted line, the dot- 
short-dashed line and the short-dashcd-long-dashed line represent 
the PDF calculated by decreasing the continuum level by 5, 3 and 
1.5 per cent, respectively. In the lower panel, and following the 
same legend as above, the various lines display the ratio between 
the PDF calculated by increasing or decreasing the continuum 
and the PDF calculated with the final continuum. The errors on 
the ratios in the lower panel are similar for all curves but plotted 
only once for clarity. 



uum placement on the PDF will be local (i.e. it should cause 
small variations around particular flux values), and will not 
affect its global average as strongly as Fig. 8 suggests. Re- 
gardless, however, when we turn to obtaining constraints 
on the thermal state of the IGM from our PDF measure- 
ments later, we follow Bolton et al. (2008) and Viel, Bolton 
& Haehnelt (2009) and only consider the PDF in the flux 
range F = [0.1-0.8]. This ignores the part of the PDF which 
is the most sensitive to continuum uncertainties, F > 0.8. 



tt 



Note this is not necessarily true for moderate to low resolution 



and S/N data, which are often normalised assuming a power- 
law extrapolation from the relatively unabsorbed spectral regions 
uncontaminated by emission lines redwards of the QSO Lyman 
a. However, more accurate approaches also for this regime, based 
on principal component analysis, have been developed (e.g. Lee 
et al. 2012, Paris et al. 2011). 
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Figure 9. The redshift evolution of the PDF with the full sample divided into two redshift bins, (z) = 2.90 (left panels) and (z) = 3.45 
(right panels). Upper panels: The squares with error bars represent the PDF measured after removing the metal line absorption. The 
dashed lines join the squares and are drawn in order to guide the eye. The solid lines are the PDFs computed for the full sample, but 
including metal absorption lines. Lower panels: The ratios between PDF calculated before and after removing the metal lines. 



3.4 Redshift evolution of the PDF 

Fig. 9 shows the PDF computed for the full sample, di- 
vided into two different redshift bins, centered at redshifts 
2.90 and 3.45. This figure also illustrates the effect of metal 
absorption lines on the PDF. The ratio between the PDF 
computed before and after removing metal absorption lines 
depends weakly on redshift. 

The PDF for the full sample, divided into two redshift 
bins, is provided in Table 2. We tabulate the PDF com- 
puted for (1) all pixels, (2) after removing contamination 
from metal absorption lines and (3) after removing both 
metal absorption lines and Lyman limit systems; the latter 
has been used for the comparison with simulation results as 
described in Sect. 4. In Fig. 10, the PDFs excluding met- 
als are compared to the PDFs obtained by K07 at lower 
redshift. For purposes of clarity, and to better focus on the 
evolution of the PDF, here the error bars are omitted. As al- 
ready noted, the spectra at higher redshift are characterised 
by lower average flux levels. In terms of the PDF, this trans- 
lates into a greater number of pixels with a low fluxes, re- 
sulting in a shallower distribution toward higher redshift. 
This trend is clearly visible in Fig. 10. 

It is also worth noting that a different method to re- 
move metal absorption lines was used by K07. In the work 
of K07, pixels belonging to spectral regions including met- 
als were not excised. The PDF was instead calculated by 
using the spectra fitted with VPFIT, including only Lyman 
a absorption lines. The use of these different methods could 
cause small differences between the PDFs. However, we have 



checked that these differences are within the error bars esti- 
mated with the jackknife method, and the results shown in 
Fig. 10 are unlikely to be significantly affected by the differ- 
ent methods used to remove metals from the QSO spectra. 

In addition, in Fig. 10 there is a difference between our 
PDF measured in the lowest redshift bin and the PDF ob- 
tained by K07 in their highest redshift bin, i.e. at (z) = 2.94. 
This difference is better illustrated in Fig 11, where we re- 
port our PDF calculated at (z) = 2.90 and the PDF mea- 
sured by K07 at (z) = 2.94. Our PDF is higher than the K07 
measurement by a factor ranging from 1.1 to 1.5 in the flux 
range 0.05 < F < 0.8. Note that two of the QSO spectra 
used in this study were also used in the PDF measurements 
presented by K07 (Q0055-269 and PKS2126-158). These ob- 
jects contribute significantly to the highest bin of K07, as 
well as to our lowest redshift bin. As a sanity check, we have 
re-calculated the PDF with the two spectra of Q0055-269 
and PKS2126-158 only, by using the pixels in the same wave- 
length range adopted by K07. The results of this test are 
visible in Fig. 12, along with jackknife error bars: once the 
same wavelength ranges as those of K07 for Q0055-269 and 
PKS2126-158 are taken into account, the two PDFs are con- 
sistent with each other within the uncertainties. This result 
tells us that the discrepancy between our PDF and K07's is 
not due to a mistake in the computation of the PDF. 

Since we expect that the flux PDF in a given redshift 
range does not depend on the sample of QSO spectra used 
to compute it, the observed discrepancy could be due to 
the modest size of the sample. If this is the case, our result 
should be more representative of the true value of the PDF 
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Figure 10. The PDF measured in this work at (z) = 2.9 (squares 
and long-dashed line), (z) = 3.45 (circles and short-dashed line) 
compared with the PDF obtained by Kim et al. (2007) at (z) = 
2.07 (light grey), (z) = 2.52 (grey) and (z) = 2.94 (dark grey). 



Figure 11. The PDF measured in this work in the lowest red- 
shift bin (squares and long-dashed curve) plotted with error bars 
compared to the PDF obtained by Kim et al. (2007) at (z) = 2.94 
(dark grey curve). 



since the sum of the Az of our QSO spectra contributing to 
the considered bin is larger than ~ 50 per cent with respect 
to K07. We note also that the average redshift of the pixels 
belonging to our lowest redshift bin is larger than the corre- 
sponding value for K07s highest redshift bin ({z p i xe i) = 3.03 
versus 2.91, respectively). This difference could contribute 
to the fact that, as seen in Fig. 10, our PDF at (z) — 2.90 is 
higher than K07s PDF at at (z) = 2.94. 

Summarizing, while the evolution with redshift of the 
flux PDF is evident from Fig. 10, more observations of QSOs 
at z ~ 3.3 could be needed to establish the true level of the 
flux PDF atz~ 2.9. 

3.5 Redshift evolution of the HI optical depth 

A quantity of interest in cosmology is the effective optical 
depth of the Lyman a forest, which is also sensitive to the 
thermal and ionization state of the IGM at high redshift. 
If (F)m = exp (— r G ff ) represents the average normalized 
flux measured from the Lyman a forest, the effective opti- 
cal depth may be defined as T e g = — ln(_F)Hi- The effective 
optical depth is therefore related to the average flux in all 
pixels in the Lyman a forest, and provides different and 
complementary measure to the PDF, which is a differential 
measurement. 

In order to better probe the redshift range investigated 
in this paper and to increase the amount of information 
which can be extracted from our sample, we have divided 
each of the eight spectra analysed in this work into two red- 
shift bins. In each bin we have then calculated the effective 
optical depth. In Tab. 3, we show the effective optical depth 
for our sample, computed both before and after removing 
metal contamination. We have calculated T e g by consider- 
ing the flux in the wavelength ranges reported in Tab. 3, 




0.2 0.4 0.6 0.8 



Flux 

Figure 12. The flux PDF measured by K07 at (z) = 2.94 
(dark grey curve) plotted with error bars compared to the PDF 
measured from the two spectra in our sample, Q0055-269 and 
PKS2126-158, used in the K07 measurement (open squares and 
long dashed curve) . This comparison uses pixels in the same wave- 
length range adopted K07. 
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Figure 13. The redshift evolution of the Lyman a forest effective 
optical depth. The solid circles are the measurements from the 
present work. The triangles are the data from Schaye et al. (2003) 
and the diamonds are from Kim et al. (2007). The dashed line 
is the best fit power-law and the solid line the power-law plus 
Gaussian which best fits the whole sample of data (see text for 
details). The dotted line is the best fit power-law plus Gaussian 
obtained by Faucher-Giguere et al. (2008) from the study of a 
sample of high-resolution spectra. 

after removing a fraction of pixels as described in section 
2.2 in order to avoid the proximity effect. Note also that 
since we wish to measure the HI effective optical depth, re- 
moving metal lines is again very important for obtaining a 
precise estimate of r e a which can be compared to other sets 
of data or with simulations. 

The errors on r e g are computed by means of a bootstrap 
method, as described in McDonald et al. (2000), Schaye et 
al. (2003) and K07. We have divided each spectrum into N 
chunks of 100 pixels, corresponding to 5 A. A total of 500 
bootstrap realizations were then performed by randomly se- 
lecting the N chunks with replacement. The standard de- 
viation from the mean r e ff is the error reported in Tab. 3. 
In the last column of Tab. 3 we also give the percentage 
contribution of metal absorption lines to T e g, which can be 
calculated as: 

metals (per cent) = ( Tcff - HI + z _ i) . 10 o. (4) 

T e ff 

The metal contamination varies from 1 to 28 per cent; we 
will discuss its redshift dependence in more detail later on 
in this section. 

In Fig. 13, we show the evolution of Lyman a forest 
optical depth, r e g, measured in this work after the removal 
of metal lines, compared to previous estimates by various 
authors. The dotted line is the best fit obtained by Faucher- 
Giguere et al. (2008) from a sample of 86 high-resolution, 
high-S/N ESI and HIRES quasar spectra. These authors 
also find evidence for a deviation of r e g from a power-law at 
z ~ 3.2, consistent with the previous result of Bernardi et 



Figure 14. Redshift evolution of the ratio between the optical 
depth computed without removing the metal absorption and af- 
ter the removal of the metal lines. The solid circles are the val- 
ues computed in the present work. The triangles are the data by 
Schaye et al. (2003) and the diamonds are the data by Kim et al. 
(2007). The dashed line represents a straight-line fit to the data, 
y = ax + 6, with a = 1.37 and b = —0.076. 

al. (2003) obtained using SDSS data and a different method- 
ology (see also Paris et al. 2011). The solid line is our best 
fit to all of the data (i.e. those of Schaye et al. 2003, Kim 
et al. 2007 and the present work) obtained by means of the 
relation 

r aB = A(l + z) B + Cexp{- [(1 + 2 ^~ ° ]2 }, (5) 

which is a modified power-law including a Gaussian com- 
ponent to take into account the presence of a deviation, as 
in Faucher-Giguere et al. (2008) . In comparison, the dashed 
line represents the best fit to the whole data ensemble for a 
power-law only. The reduced \ 2 are reported in Fig. 13; we 
find that a power law provides a better fit to the data. In the 
case of a power law plus Gaussian, the best fit was obtained 
with A = 0.0023, B = 3.64, C = -0.05, D = 4.36, E = 0.21, 
i.e. with parameters associated with a very tiny deviation 
from a power law. Note also that the reduced \ 2 reported in 
Fig 13 are in general larger than the values reported in table 
5 of Faucher-Giguere et al. (2008). However, the data used 
by Faucher-Giguere et al. (2008) are more homogeneous and 
show a lower dispersion, hence it seems reasonable for their 
fits to show lower \ 2 values. From these results, we conclude 
that, from our study of the redshift evolution of the effec- 
tive optical depth, we do not find any strong evidence for 
a deviation from a power law. The best fit is expressed by 
the power law T cff = 0.0028 ± 0.0003 (1 + z ) 3 45±0 08 j i n ver y 
good agreement with the previous best fit provided by K07: 
r cff = 0.0023 ± 0.0007(1 + z f- 65±0 - 21 . 

In Fig. 14, we show the redshift evolution of the ratio 
between the effective optical depth measured before and af- 
ter removing metal absorption lines. This plot illustrates the 
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effect of metal absorption on the Lyman a forest effective 
optical depth as a function of redshift. Our data are plot- 
ted together with the measurements from K07 and Schaye 
et al. (2003). The dashed line represents a straight-line fit 
to the data, y = ax + b. Our best fit values are a = 1.36 
and b = —0.07. The weak correlation shown in Fig. 14 again 
highlights the previously discussed redshift dependence of 
the metal line contamination in the Lyman a forest, which 
is consistent with a greater amount of metal absorption to- 
ward lower redshift. 



4 THE SIMULATIONS 

In this section, we now turn to describe the suite of hydrody- 
namic simulations used to obtain constraints on the thermal 
state of the IGM at redshift 0^3 from our PDF measure- 
ments. Although the basic setup of the simulations is similar 
to earlier work by Bolton et al. (2008) and Viol et al. (2009), 
we now use new, higher resolution runs to explore a wider 
thermal history parameter space compared with those used 
by Viel et al. (2009). The simulations were performed us- 
ing the parallel Tree-PM smoothed particle hydrodynamics 
(SPH) GADGET-2 code (Springel 2005) and its recently up- 
dated version GADGET-3. All the simulations were started 
at z = 99, with initial conditions generated using the trans- 
fer function of Eisenstein & Hu (1999). The gravitational 
softening length was chosen to be equal to l/30th of the 
mean linear interparticle spacing and star formation was 
included using a simplified prescription which converts all 
gas particles with temperature T < 10 5 K and overdensity 
A = p/{p) > 10 3 into collisionless stars. 

In Table 4 we summarize the parameters used in the 
cosmological simulations. We use two different sets of sim- 
ulations. Set A includes high resolution runs and, among 
them, the fiducial run for this paper: AIref- Each simula- 
tion has a box size of 10 h~ Mpc comoving and contains 
2 x 512 3 dark matter and gas particles. The model chosen 
is a flat ACDM cosmology with parameters: ilom = 0.26, 
n A = 0.74, Q b h 2 = 0.023, n s = 0.96, H = 72 km s" 1 
Mpc -1 and as = 0.80, which are in agreement with recent 
studies of the cosmic microwave background (Komatsu et 
al. 2009; Reichardt et al. 2009; Jarosik et al. 2011). 

We use set A to explore variations in the thermal 
state of the gas. A spatially uniform ultra-violet background 
(UVB) produced by quasars and galaxies (Haardt & Madau 
2001) was adopted, where hydrogen is reionized at 2 = 9. 

The gas in the simulations is assumed to be optically 
thin and in ionization equilibrium with the UVB (i.e. photo- 
heating by the UVB is balanced by cooling due to adiabatic 
expansion). In these physical conditions, gas with overden- 
sity A ^ 10 follows a tight power-law temperature-density 
relation: T = To A 7-1 ( Hui & Gnedin 1997; Valageas, 
Schaeffer & Silk 2002), where To is the temperature of the 
IGM at the mean density. In order to obtain various thermal 
histories, the Haardt & Madau (2001) photo-heating rates 
eP M01 were rescaled by different values: e; = kA q e? M01 , 
where i =[HI, Hel, Hell]. Our fiducial simulation AIref 
has parameters k = 2.20 and q — —1.00 (which produces an 
isothermal T — A relation). Runs AIcld and AIhot also 
have q = —1.00 (still an isothermal T — A relation), but 
k — 1.45 (producing lower To) and k = 3.10 (higher To), 



respectively. Runs AIlw-y an d Alm-y both have k — 2.20, 
but q = —1.60 (resulting in an inverted T — A relation with 
7 = 0.7) and q = -0.45 (flattened T - A relation with 
7 = 1.3), respectively. 

Simulation set B includes lower resolution runs that are 
intended to explore the dependency of the PDF on cosmo- 
logical parameters. These simulations have a box size of 10 
h^ 1 Mpc comoving and contain 2 x 256 3 dark matter and 
gas particles. Again, the reference cosmological model (see 
run BIref in Tab. 4) is a flat ACDM model, but now the 
parameters are slightly different: Qom = 0.26, Qa = 0.74, 
n b h 2 = 0.024, n s = 0.95, H = 72 km s _1 Mpc -1 and 
as — 0.85. As shown in Tab. 4, the cosmological parameters 
varied are as (runs BIlss and BIhss), n s (runs BIlns and 
BIhns), fJom (runs Bl L Om and Bl H Om) and H (runs BIlho 
and BIhho)- Note that this second set of simulations has a 
slightly different as value than the first set, however, we do 
take this into account in the analysis described below. We 
refer the reader to Becker et al. (2011), Bolton et al. (2008) 
and Viel, Bolton & Haehnelt (2009) for further details on 
the simulations. 



5 ANALYSIS OF THE SIMULATIONS 
5.1 Synthetic spectra and PDFs 

Given the positions, velocities, densities and temperatures of 
all the SPH particles at a given redshift, spectra along lines- 
of-sight (LOSs) through the simulation boxes were com- 
puted following the procedure described by Theuns et al. 
(1998). The interested reader can find more details about 
this procedure in Section 5 of Tescari et al. (2011). Two dif- 
ferent sets of synthetic spectra were constructed: the first 
set includes two redshift bins at z = 2.95 and 3.48 and 
the second set includes three redshift bins at z = 2.90, 3.25 
and 3.55. After extracting the spectra along random LOS 
through the cosmological box at redshift z, we rescaled all 
the HI optical depths by a constant factor, so that their 
mean value was equal to the HI effective optical depth, r e s, 
given by the Kim et al. (2007) fit: 

r cB = (0.0023 ± 0.0007) (1 +z) 3 - 65±0 ' 21 . (6) 

This rescaling ensures that our spectra match the observed 
mean normalized flux of the Lyman a forest at the appro- 
priate redshift: (-F , ) H i,obs = exp (— T e s). The spectra are then 
convolved with a Gaussian of 7 km s _1 FWHM and rebinned 
on to pixels of width 0.05 A. Finally, in order to have realis- 
tic spectra to compare with observations, we add Gaussian 
distributed noise with total signal to noise/pixel (i.e. stan- 
dard deviation of the Gaussian at F — 1) equal to 60 and 
read out S/N equal to 100. 

The dependence of the simulated flux PDF on To, 7 and 
as is illustrated in Fig. 15. The two upper panels show the 
PDFs extracted from our two fiducial simulations: the high 
resolution run AIref (black solid curve) and the low reso- 
lution run BIref (red dashed curve) at redshift z — 2.95 
(left-hand side) and z — 3.48 (right-hand side). Each model 
is compared to the PDF of the full observational sample 
measured after the removal of metals and LLS (data points 
with error bars). The simulated PDFs in this instance have 
been obtained by rescaling the value of the effective op- 
tical depth, r ff, in order to minimize the \ 2 statistic. In 
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Table 4. Summary of the hydrodynamical simulations used in this work. Column 1, simulation name: note the distinction between 
subset A (high resolution set) and the subset B (low resolution set); column 2, number of gas particles; column 3, mass of each gas 
particle; column 4, matter density parameter. In all models we assume a flat universe with = 1 — S7om- However, Q^h 2 = 0.023 for 
the subset A and U h h 2 = 0.024 for the subset B of our simulations; column 5, primordial power spectrum index; column 6, Hubble 
constant; column 7, fluctuation amplitude at a scale of 8 h^ 1 Mpc; column 8, slope of the temperature-density (T — A) relation, where 
T(A) = To A 7_1 ; column 9: gas temperature at the mean density, in units of 10 3 K, at z = 3.25; column 10: further details on the 
simulations. 



Run 




m-GAS 






H 




7 


T (10 3 K) 


Notes 






(h- 1 ) M 






(km/s/Mpc) 






(z = 3.25) 




AIref 


512 3 


9.2 x 10 4 


0.26 


0.96 


72 


0.80 


1.0 


19.42 


Fiducial high res. model 




512 3 


9.2 x 10 4 


0.26 


0.96 


72 


0.80 


1.0 


14.93 


(~!nlH run 


AIhot 


512 3 


9.2 x 10 4 


0.26 


0.96 


72 


0.80 


1.0 


23.88 


Hot run 


AI.LW7 


512 3 


9.2 x 10 4 


0.26 


0.96 


72 


0.80 


0.7 


19.96 


Low 7 


A1hi 7 


512 3 


9.2 x 10 4 


0.26 


0.96 


72 


0.80 


1.3 


19.07 


High 7 


BIref 


256 3 


7.4 x 10 5 


0.26 


0.95 


72 


0.85 


1.0 


23.79 


Low res. reference model 


BIlss 


256 3 


7.4 x 10 5 


0.26 


0.95 


72 


0.80 


1.0 


23.64 


Low as 


B1hS8 


256 3 


7.4 x 10 5 


0.26 


0.95 


72 


0.90 


1.0 


23.95 


High o- 8 


BIlns 


256 3 


7.4 x 10 5 


0.26 


0.90 


72 


0.85 


1.0 


23.67 


Low n B 


BIhns 


256 3 


7.4 x 10 B 


0.26 


1.00 


72 


0.85 


1.0 


23.93 


High n s 


BlLOm 


256 3 


7.4 x 10 5 


0.22 


0.95 


72 


0.85 


1.0 


24.55 


Low f2om 


BlHOm 


256 3 


7.4 x 10 5 


0.30 


0.95 


72 


0.85 


1.0 


23.12 


High n 0ul 


BIlho 


256 3 


7.4 x 10 5 


0.26 


0.95 


64 


0.85 


1.0 


25.40 


Low Ho 


BIhho 


256 3 


7.4 x 10 5 


0.26 


0.95 


80 


0.85 


1.0 


22.35 


High H 



the analysis presented in Sec. 5.2, we will only compare our 
simulations to flux values in the range [0.1 — 0.8], which is 
less sensitive to possible continuum errors. To demonstrate 
this, the simulated PDFs from run AIref (blue dot-dashed 
lines) and BIref (orange triple dot-dashed lines) with ad- 
ditional ±1.5% (left-hand side) and ±3% (right-hand side) 
continuum errors are overplotted. Even in the worst case of 
±3% errors, the simulated PDFs are consistent in the range 
F = [0.1 -0.8]. 

The lower six panels of Fig. 15 show the ratio be- 
tween PDFs obtained from simulations with different val- 
ues of astrophysical and cosmological parameters, and the 
corresponding reference PDFref (calculated by using either 
AIref or BIref), at the same two redshifts. We vary the 
following parameters: To (first row), 7 (second row) and as 
(third row). It is clear from the figure that variations of 7 
have the largest effect on the simulated PDF (see also Fig. 
2 of Bolton et al. 2008 and associated discussion). 



5.2 Constraints on cosmological and astrophysical 
parameters 

We next use the synthetic spectra to construct flux PDFs 
which may be compared to our observational measurements. 
We explore the effect of varying the following cosmological 
and astrophysical parameters on the PDF: as, n B , Qom and 
Ho for the cosmological parameters, and To and 7 for the 
IGM thermal history. Note, however, that most of the cos- 
mological parameters are only weakly constrained by the 
flux PDF. The largest effect is induced by using a differ- 
ent value of as (e.g. Bolton et al. 2008). We will therefore 
present marginalized cosmological constraints for the ampli- 
tude of the matter power spectrum only. The astrophysical 
parameters To and 7 are parameterized as a broken power- 
law around z = 3, in order to investigate possible changes 
in the thermal state around Hell reionization: To = T^(z = 



3) x [(1 + z)/4] T o« and 7 = -y A (z = 3) x [(1 + z)/4p S M 
(Viel et al. 2009). The effective optical depth evolution is 
parameterized r cff = t^(z = 3) x [(1 + z)/4] T rff (z) . 

In order to explore our parameter space efficiently with 
our limited number of simulations, we Taylor expand the 
flux PDF and compute derivatives to second order using 
two simulations (plus the reference run) for each parame- 
ter. This closely follows the procedure described in Viel & 
Haehnelt (2006), where this Taylor expansion method was 
used in order to explore constraints derived from the Sloan 
Digital Sky Survey flux power spectrum (McDonald et al. 
2006). We start with the AIref simulation and proceed to 
calculate the \ 2 °f our models by varying all the parame- 
ters that alter the flux PDF, enabling us to expand around 
the fiducial starting model. The expansion method allows us 
to explore the parameter space close to the reference model 
with an accurate set of hydrodynamic simulations, even if 
the full parameter space cannot be probed in this way with 
the same high accuracy. We impose weak priors on the ef- 
fective optical depth evolution, but both its amplitude and 
slope are constrained strongly by the data, so the priors 
do not affect our results. In order to present a conservative 
analysis, we also decide to consider only the values of the 
flux in the range F = [0.1 — 0.8], as this range should not 
be affected much by the errors made in the continuum fit- 
ting procedure (e.g. K07; Bolton et al. 2008, and see Fig. 8). 
Furthermore, when estimating the covariance between the 
PDF data points we use the covariance extracted from the 
simulations for the non-diagonal terms, which is less noisy 
(e.g. Lidz et al. 2006). 

The one dimensional marginalized constraints are 
shown in Figure 16: the likelihoods obtained in the present 
work are represented by the black curves, while the results 
derived from z < 3 K07 data set analysed by Viel et al. 
(2009) are shown as dashed blue lines. The joint constraints 
for both PDF measurements are plotted as the red lines. 
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Figure 15. Simulated PDFs at redshift z = 2.95 (left-hand side) and z = 3.48 (right-hand side). Upper two panels: Our two reference 
simulations are displayed in both plots: run AIref (black solid curves) and run BIref (red dashed curves). The black data points 
show the PDF measured from the full observational sample after the removal of metals and LLS (Tab. 2, last two columns). The blue 
dot-dashed curves and orange triple dot-dashed curves are the PDFs obtained with ±1.5% (left-hand side) and ±3% (right-hand side) 
continuum errors added to the spectra from the runs AIref and BIreFj respectively. The two vertical dotted lines mark the region inside 
which we perform the fit to the observational data in this work. Lower six panels: the ratio between PDFs obtained from simulations 
with different values of astrophysical and cosmological parameters, and the corresponding reference PDFref- Variations in the following 
parameters are explored: To (first row), 7 (second row) and ag (third row). 



Table 5. 1 — a parameter constraints from the sample of QSOs analysed in a) this work and b) the joint sample composed of this work 
and the K07 PDF data. 





T A 

off 


T S 
off 


7 A 7 s (z < 3) 


7 s (2 > 3) 


T A [K] 


T$ (z < 3) 


T$ (z > 3) 




this work 


0.322 


3.15 


0.90 


-1.8 


19250 




-0.8 


0.8 




±0.012 


±0.28 


±0.21 


±1.0 


±4800 




±1.5 


±0.04 


this work - 


- K07 data 0.308 


3.26 


0.70 -1.4 


-0.6 


17900 


0.1 


-0.1 






±0.005 


±0.10 


±0.12 ±0.7 


±1.3 


±3500 


±1.0 


±1.6 
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Figure 16. One-dimensional marginalized constraints on the cosmological and astrophysical quantities used in the analysis (see text for 
details of the parametrization). The blue curves indicate the constraints obtained by using the PDF measurements at z < 3 presented by 
K07 (see also the analysis of Viel, Bolton & Haehnelt 2009), black-continuous curves refer to the present work, while the joint constraints 
for both this data set and the measurements of K07 are shown as continuous red curves. 



From this figure we can see that overall there is agreement 
between the inferred parameter values for the high redshift 
(this work) and low redshift (K07) measurements. The 1 — a 
parameter constraints from the sample of QSOs analysed in 
this work and for the joint sample including the K07 data are 
summarized in Tab. 5. The minimum \ 2 f° r the best fitting 
model to the PDF measured in this work is 52 for 38 degrees 
of freedom (d.o.f.) which is likely to happen 7 per cent of 
the time. In contrast, for the joint sample the minimum \ 2 
is 94 for 83 d.o.f. with has a 21 per cent probability. The 
probability is calculated using the full covariance matrix for 
all the parameters recovered. Note that the "degrees of free- 
dom" value does not have a precise meaning if the error bars 
are heavily correlated. However, this happens not to be the 
case in the flux ranges considered here, F = [0.1-0.8], where 
correlations are weak. We have also checked that there is 
good agreement between the derived constraints from the 



present data set when compared to constraints obtained by 
splitting the sample in two and three redshift bins. 

5.3 Implications for the IGM thermal history 

Finally, in Figure 17 we display the redshift evolution of 
the temperature at the mean density, To, and the slope of 
the temperature-density relation, 7, inferred from our anal- 
ysis. We use the broken power-law parameterization around 
z = 3 introduced previously: / = f A x [(l + z)/4]^ where 
/ = To or 7, / A is the corresponding normalization factor at 
z = 3 and f S (z) is the slope of the relation (see Tab. 5). The 
black solid lines in Figure 17 are the results obtained from 
this work (z > 3) while the shaded blue regions show their 
associated errors. The red dot-dashed lines are the results 
obtained from the joint sample composed of this work and 
the K07 data. The red dashed lines show the associated er- 
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Figure 17. Redshift evolution of the temperature at the mean 
density, To, (upper panel), and the slope of the temperature- 
density relation, 7, (lower panel) inferred from our analysis of 
the flux PDF. The black solid lines are the results obtained from 
this work (z > 3), with the shaded blue regions showing the asso- 
ciated errors (see text for details). The red dot-dashed and lines 
are the results obtained considering the joint sample composed 
of this work and the K07 data. The associated errors are shown 
by the red dashed lines. The vertical black dotted lines mark the 
transition between z < 3 and z > 3. 



curvature statistic (Becker et al. 2011) are typically better 
suited for this measurement. 

On the other hand, the flux PDF is sensitive to the slope 
of the temperature-density relation, 7. The constraints dis- 
played in Fig. 17, based on the PDF measurements presented 
in this work, are formally consistent with an isothermal 
temperature-density relation at z = 3 (i.e. 7 ~ 1), although 
the uncertainties remain large. A temperature-density rela- 
tion which is close to isothermal may be expected if Hell 
reionization was underway by z = 3 (e.g. Schaye et al. 2000, 
Theuns et al. 2002). When combining our data with the K07 
measurements at lower redshift, however, the inferred 7 at 
2 = 3 is shifted to lower, inverted (7 < 1) values, consistent 
with previous results based on the K07 data (e.g. Bolton 
et al. 2008, Viel et al. 2009). This is because the K07 data 
is more constraining at z < 3 due to the larger number 
of QSOs in this sample, and for the simple power-law pa- 
rameterisation we have assumed this also pushes our joint 
constraints to smaller 7 values at z > 3. 

It has been suggested recently that an inverted 
temperature-density relation at may be explained by an 
IGM which is volumetrically heated by TeV emission from 
blazars (e.g. Chang et al. 2011, Puchwein et al 2011). This is 
in contrast to conventional photo-ionization heating models 
which predict a temperature-density relation which evolves 
from isothermal to 7 ~ 1.6 following reionization (Hui & 
Gnedin 1997). We caution, however, the uncertainties on 
the flux PDF measurements are still large; the temperature- 
density relation suggested by the radiative transfer simula- 
tions of McQuinn et al. (2011), 7 ~ 1.3 at z = 3, is com- 
patible at the 1.5(7 level with the value derived from the 
sample considered in this work alone. The K07 data alone 
are furthermore consistent with an isothermal IGM, again to 
within 1.5a (Viel et al. 2009). We conclude that tighter con- 
straints on 7 from larger data sets, as well as a comparison of 
different measurement techniques to aid with the identifica- 
tion possible systematics, will be required for distinguishing 
between different IGM heating scenarios. 



rors. Since the / A and f s (z) parameters are correlated and 
the f S (z) are in general poorly constrained, we define the 
error associated to /, at a given redshift z, as follows: 

A/=(/ u V^)x(l±£) /SW , (7) 

where / A = f A + la and /jj = / A — la are, respectively, 
the 1 — a upper and lower limits on / A , and f s (z) are the 
best fit values of the slope showed in Tab. 5. 

The temperature at mean density and its redshift evolu- 
tion, shown in the upper panel of Figure 17, are not strongly 
constrained by the flux PDF, although our measurements 
are consistent with other recent analyses presented in the 
literature which hint that Hell reionization may be com- 
pleting around z ~ 3 (Lidz et al. 2010, Becker et al. 2011). 
However, the PDF constraints alone are not precise enough 
to draw conclusions about the timing and extent of Hell 
reionization; the flux PDF is a relatively insensitive to To 
(see e.g Fig. 15). Alternative measurements from Doppler 
widths (Schaye et al 2000) wavelets (Lidz et al. 2010) or the 



6 CONCLUSIONS 

In this paper, we have measured the Lyman 01 forest flux 
PDF from a set of eight QSO spectra to investigate the 
thermal state of the IGM at z ~ 3. The absorption lines 
in the Lyman a forest have been fitted with Voigt profiles, 
and absorption lines arising from heavy elements have been 
carefully identified. The spectra have then subsequently had 
the metal absorption lines removed from within the Lyman 
a forest. The resulting masked spectra have been used to 
study the redshift evolution of the Lyman a forest flux dis- 
tribution and the evolution of the HI effective optical depth. 
Our results can be summarized as follows. 

• By removing metals from the Lyman a forest, we at- 
tempt to eliminate pixels which will introduce a systematic 
uncertainty into our measurement of the flux PDF. The ef- 
fect of removing these metals on the PDF is to decrease the 
number of pixels in bins of low flux and increase the number 
of pixels in bins of high flux. As a result, the PDF measured 
after correcting for metal absorption is slightly steeper than 
the PDF measured without removing metals. In general, the 
effect of metal removal on the PDF is strongest for flux bins 
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at F < 0.6, in agreement with previous results for other 
QSO samples at lower redshift (K07). 

• We measure the PDF for our full sample in two redshift 
bins, both before and after metal decontamination, and find 
the effect of metal absorption on the Lyman a forest flux 
PDF is redshift dependent. It is more significant at lower 
redshift, where the average opacity of the Lyman a forest is 
lower and there are fewer absorption lines. This redshift de- 
pendence may be understood in terms of the larger number 
of metal lines which are blended with Lyman a absorption 
lines toward higher redshift, resulting in the identification 
of fewer isolated metal absorption lines. On the other hand, 
it is also possible that the metallicity of the intergalactic 
medium decreases with increasing lookback time. 

• We present measurements of the effective optical depth, 
T e ff, from our QSO sample, which provides a probe of the 
IGM which is complementary to the PDF. The values of r e g 
measured from the QSOs in our sample can be fitted with 
a power-law r cff = 0.0028 ± 0.0003 (1 + z ) 3 - 45±0 08 ; in agree- 
ment with previous measurements by K07 after correcting 
for metal absorption. We also find no substantial evidence 
for a deviation from a power law in the log(r e ff ) — log(l + z) 
plane at z ~ 3.2. 

There is a weak redshift dependence of the ratio of the HI 
effective optical depth measured before and after removing 
identified metal absorption lines. The ratio increases toward 
lower redshift, and is consistent with our previous finding 
from the PDF that metal absorption plays a more signifi- 
cant role in the Lyman a forest toward lower redshift. 

• We perform a set of high resolution hydrodynamical 
simulations with a range of cosmological and astrophysical 
parameters and use them to extract mock QSO spectra with 
properties that closely resemble those of the observed spec- 
tra. We compute the Lyman a flux PDF and compare it to 
the observed PDF in the range F = [0.1-0.8] via a Taylor 
expansion at second order to find deviations from a best- 
guess case (e.g. Viel & Haehnelt 2006). By implementing 
this method within a Monte Carlo Markov Chain likelihood 
estimator, we present constraints for the IGM thermal state. 
At 2 = 3, our PDF measurements are consistent with an 
isothermal temperature-density relation, 7 = 0.90 ± 0.21, 
and a temperature at the mean density of To = 19250 ±4800 
K (la uncertainties). 

• We also obtain joint constraints on the IGM thermal 
state by analysing the lower redshift PDF measurements 
of K07, which is a significantly larger data set than the 
one considered here. At z = 3, we find the joint data set 
is consistent with an inverted temperature-density relation, 
7 = 0.70 ± 0.12, with To = 17900 ± 3500 K (la uncertain- 
ties). 

• In conclusion, the present work highlights the role of 
the flux PDF as a competitive and quantitative thermome- 
ter which may be used to constraint the IGM temperature- 
density relation over a wide redshift range. The constraints 
derived from the present work are complementary with re- 
spect to other approaches, such as the wavelets method (Lidz 
et al. 2010), the cutoff in the Doppler width-column density 



plane (Schaye et al. 2000) or the curvature (Becker et al. 
2011), which are primarily sensitive to higher density re- 
gions of the IGM at z ~ 3. In the near future, the most 
promising route to making further progress on our knowl- 
edge of the thermal state of the IGM will be to carefully 
compare and contrast results from these different probes of 
the IGM temperature-density relation. 
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